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We study the late stages of spinodal decomposition in a Ginzburg-Landau mean field model with 
quenched disorder. Random spatial dependence in the coupling constants is introduced to model 
the quenched disorder. The effect of the disorder on the scaling of the structure factor and on the 
domain growth is investigated in both the zero temperature limit and at finite temperature. In 
particular, we find that at zero temperature the domain size, R(t), scales with the amplitude, A, of 
the quenched disorder as R(t) — A~@ f(t/ A -1 ) with (5 ~ f .0 and 7 ~ 3.0 in two dimensions. We 
show that /3/7 = a, where a is the Lifshitz-Slyosov exponent. At finite temperature, this simple 
scaling is not observed and we suggest that the scaling also depends on temperature and A. We 
discuss these results in the context of Monte Carlo and cell dynamical models for phase separation 
in systems with quenched disorder, and propose that in a Monte Carlo simulation the concentration 
of impurities, c, is related to A by A ~ c 1//d . 

PACS: 64.60.-i, 64.60.Ak 



I. INTRODUCTION 

Spinodal decomposition has received a great deal of 
attention in recent years since it is relevant to many ma- 
terials processing problems and is one of the unsolved 
problems in statistical mechanics. Much progress has oc- 
curred over the past decade in understanding both the 
early and late time behavior of spinodal decomposition 
in the simplest systems with no disorder and a single 
order parameter but there is more to learn con- 

cerning the important process of spinodal decomposition 
in disordered systems. Systems with quenched disorder, 
for example, are interesting because the kinetics of spin- 
odal decomposition differs markedly from the behavior of 
pure systems. In particular, the usual power law growth 
of the characteristic domain size, R ~ t a is thought to 
change over to a logarithmic growth Q 
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The earliest computational studies of phase separation 
in systems with quenched disorder were Monte Carlo 
[MC] simulations performed nearly a decade ago |^| ||. 
MC simulations are generally inefficient for studying the 
domain growth since diffusion of particles through the 
system at late stages requires many updates. However, 
MC simulations have the distinct advantage of introduc- 
ing the quenched disorder in a very natural way. Using 
MC simulations, Crest and Srolovitz |(| found evidence 
of a scaling function relating R to the concentration of 
impurities in a two-dimensional Ising model with non- 
conserved order parameter, but were unable to confirm 
the prediction Q of logarithmically slow domain growth, 
Eq. ( |l.ip ■ Similar behavior was seen in MC simulations of 
a model for simultaneous phase separation and gelation 
in an Ising-type model with conserved order parameter 



Further progress has recently been made ||10| , pd| | us- 
ing more efficient cell dynamical methods |12[|. These 
cell dynamical simulations were able to probe the late 
stage behavior for both cases of conserved and noncon- 
served order parameter by introducing the disorder as 
a phcnomcnological parameter in the Ginzburg-Landau 
free energy Using cell dynamical methods, Puri ct 

ai. 0,y| were able to observe the logarithmically slow 
domain growth Q for both conserved and nonconserved 
order parameter, although the exponent in the pre- 
dicted 



behavior, 
order parameter 



Eq. (1.1), for the case of conserved 
, appeared to vary with the amplitude 
of the quenched disorder. This could be due in part to 
the systems not having reached their asymptotic behav- 
ior; even using the cell dynamical approach, simulation 
of disordered systems is computationally demanding due 
to the logarithmically slow growth. 

Here we study a time-dependent Ginzburg-Landau 
(TDGL) model and consider scaling of the domain size as 
a function of the quenched disorder both in the zero tem- 
perature limit and for T > 0. For T = 0, we expect that 
domain growth will effectively stop as various parts of the 
system fall into metastable states from which they cannot 
escape. Although the T = limit is clearly unphysical, 
domain growth pinning is seen in experimental systems 
with quenched disorder fLf] , Studying the TDGL model 
in this limit gives us the opportunity to study the effect 
of a single parameter, the quenched disorder A, on the 
dynamics of domain growth. We find scaling that relates 
the characteristic domain size R to the phenomenologi- 
cal parameter A that describes the disorder in the model. 
We also examine the behavior of this model at finite tem- 
perature. We find that this simple scaling doe s no t hold 
and propose a form for the exponent 9 in Eq. (1. 
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The rest of this paper is organized as follows. In Sec- 
tion [n], we describe the time dependent Ginzburg-Landau 
model with quenched disorder and briefly discuss its nu- 
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merical implementation. In Section III, we describe the 



results of our simulations and propose a scaling form that 
is consistent with these results. Finally in Section IV, we 



place these results in context with results for the Monte 
Carlo and cell dynamical models described above. 



II. THE MODEL AND ITS SCALING 
PROPERTIES 

In this section we briefly review the time dependent 
Ginzburg-Landau (TDGL) equation and the rescaling of 
the variables that leads to the Cahn-Hilliard equation 
|L5| , the basis of our simulation model. In addition, we 
introduce the quenched disorder into the Cahn-Hilliard 
equation [^0|JTl 1 via additional random, static coupling 
constants taken from a uniform distribution of width A. 
We present a scaling hypothesis for the model based on 
observations made in previous work JlO| , pl| . 

Consider a system with two independent species. The 
order parameter, <f>, is defined JlJ to be the difference in 
the local concentration of these two species 



: ci - c 2 . 



(2.1) 



For a conserved order parameter, the local concentration 
difference obeys a continuity equation of the form 



M <f>(r,t) - 
Here the current is 



-V-j(r,t)+C(r,f). 



j(r,t) = -MV/i(r) 



(2.2) 



(2.3) 



where M is the mobility (assumed to be constant), /i(r) 
is the chemical potential (the variational derivative of the 
free energy), and ((r, t) is Gaussian-distributed random 
noise that satisfies 

(((r,i),((r',f')) = -a B rMV 2 §{v-v) S(t-t), 

(2.4) 

where the brackets represent a thermal average over 
many different configurations. This term guarantees that 
the proper equilibrium state is achieved. In the TDGL 
model one usually chooses the Ginzburg-Landau free en- 
ergy functional 



F = J dv{ l -K\V^\ 2 -\ f 



b 2 + ~u<f> 4 } (2.5) 



Combining the c onti nuity equation, Eq. (2.2), with the 
free energy, Eq. (2^), yields 



Model B is a time-dependent Ginzburg-Landau model or 
relaxation model since the time evolution of the order 
parameter depends onl y on the minimization of the free 
energy functional, Eq. (2.5). 

The order parameter <j) is redefined by scaling t he p hc- 
nomenological parameters M,n,p and u in Eq. (2J3) as 
well as the dynamical variables r and t to obtain a simpli- 
fied equation of motion for the rescaled order parameter. 
This is known as the Cahn-Hilliard equation [ pi , and 
has the form 

|:^ = iv 2 (-V + V 3 -VV) + ^e (2.7) 

The transformation that rescales Model B into the Cahn- 
Hilliard model, following Grant et al. jl(| and Rogers et 
al. [0, is 



1/2 



r. 
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(2.8b) 



(2.8c) 



(2.8d) 



Hence from Eq. (12 



(£(x,tU(x,t)) = 



- 2k B Tu/p 2 V2 M V 2 S(x - x) 5(t - r'). (2.9) 

The formulation of model B in these variables clearly 
shows that this model involves only a single parameter, 
e, proportional to the temperature. 

To incorporate quenched disorder into Model B, it is 
conventional to introduce a spatial dependence in the pa- 
rameters p and u of Eq. (2.5) Consider couplings 
of the form p = po + 5p(x.) and u = uq + Su(x), where po 
and uq are identified with the free energy in the absence 
of disorder. After pe rfor ming the transformation of vari- 
ables leading to Eq. (2/7), the time dependent Ginzburg- 
Landau equation with quenched disorder becomes 



d_ 

dt 



1> 



dt 



MV 2 {- 



-/tV 2 c 



(2.6) 



The continuity equation, Eq. (2.2), guarantees conserva- 
tion of the order parameter <j>; this is known as Model B. 



Po 



V- 1 + 



<5it(x) 
u 



V 2 ip 



(2.10) 
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We choose for simplicity the factors 111 



Po 



*U(X) 3 



M 



(2.11) 



to be random variables taken from a uniform distribution 
of width A, thereby giving a random spatial dependence 
to the rescal ed co uplings in the model. 

Equation ( |2.10D is a generalization of the Cahn-Hilliard 
equation, Eq. ( |2.7| ). Earlier results clearly indicate that 
the crossover time, t x , from the early stage algebraic 
growth to logarithmically slow growth at later times in- 
creases as the amplitude A of the quenched disorder de- 
creases 11 . Thus it is appropriate to consider two dis- 
tinct time domains separated by a crossover region char- 
acterized by a time t x ■ For t <C t x we expect the domain 
size to grow as 



(2.12) 



where a = 1/3 from the Lifshitz-Slyosov theory p8| . 
This exponent has been observed in a wide variety of 
TDGL and cell dynamical simulations ||]4|] with con- 
served order parameter. For t t x we expect that the 
domain size at any time in the slower growth region will 
scale as R(t) ~ A~@ . Finally using the observation by 
Puri and Parekh [|ll] that the crossover time f x increases 
for decreasing disorder amplitude suggests that 



tx ~ A-< 
and hence the scaling form 

R{t)~A-f>f(t/t x )~A-Pf(tA'r), 



(2.13) 
(2.14) 
(2.15) 



From ( |2 . 1 2| ) and ( J2.14Q it follows that 

(3 = cry. 

III. CALCULATIONS AND RESULTS 



We now describe the results of our com puter simula- 
tions that verify the scaling hypothesis ( [2.14 ). First 
we calculate the average domain size, R{t), of the sys- 
tem. Since there is only one relevant length scale in a 
system undergoing spinodal decomposition W , any char- 
acteristic length, R, is proportional to the domain size. 
We choose R to be the inverse of the first moment k\ of 
the spherically-averaged structure factor, S(k,t), which 
is the structure factor S(k, t) = {J G?re _lk '( r ~ r >g(r — r')) 
averaged over the two-dimensional density of states. As 
usual, S(k, t) is the Fourier transform of the density- 
density correlation function g(r—r'). The brackets repre- 
sent the averaging over many realizations of the disorder 
parameterized by A (see Table [|) . Since we are consider- 
ing systems with conserved order parameter, the domain 
size scales as R(t) ~ l/k max , where k max is the value at 



which S(k, t) is maximum. In our simulations, however, 
we use the first moment k\ of S(k, t) (which is equivalent 
[GJ), since the peak maximum is difficult to determine 
accurately. 



We integrate the equation of motion Eq. (2.1C) directly, 
using the central difference approximation for the spatial 
derivatives. Following Rogers et al. ]l7| , we use a time 
step of At = 0.15 an d latt ice spacing Ax = 0.30 to insure 
the stability of Eq. Q2.10Q . We allow for 3 x 10 6 updates 
for each run on a two dimensional lattice of size 256 x 256 
and calculate S(k,t) at logarithmically spaced time inter- 
vals. We obtain the results discussed in this section by 
averaging over 15 to 25 realizations of A (see Table |), 
ranging from 0.0 to 0.45. 

We start with a homogeneous two-dimensional system 
with the rescaled order parameter ip initialized from a 
uniform random distribution of range —0.005 to 0.005. 
For T > 0, we use a noise amplitude y/e = 0.5, see 
Eq. ( 2.10|) . In Fig. |] we display the value of tp for every 
choice of A. Each snapshot is made after 3 x 10 6 updates 
in the time evolution. As can be seen from these snap- 
shots, the larger the disorder parameter A, the smaller 
the asymptotic domain size. The domains gradually de- 
crease in size as we increase A. This holds for simulations 
at both T = 0, the top row, and T > 0, the bottom row. 
Note that the domains have rougher and wider interfaces 
for T > 0. 

Following the discussion in Sec. 0, we calculate the do- 
main size R(t) from the first moment of the spherically- 
averaged structure factor. In Fig. ||a, we display our 
results for S(k,t) vs. wave number k = | k |, aver- 
aged over 20 runs with no quenched disorder. As ex- 
pected, p|,p7|, S(k,t) displays a peak that increases with 
time to smaller values of wave number k. The peak of the 
structure factor is larger and occurs at a greater value 
of k compared to the case with quenched disorder (see 
Fig. ||b). In fact, the finite value of the asymptotic do- 
main size for A = only exists due to the finite size 
of the lattice. This behavior is contrasted with Fig. ||d 
which shows S(k, t) for A = 0.45. As Fig. |^b suggests, 
the peak of S(k,t) reaches a maximum value at a value 
of k > 2tt/L. 

Fig. H shows the domain size R(t) at (a) T = 
and (t>) T > for each value of A. As can be seen in 
Fig. [|a, R(t) becomes constant after 10 5 — 10 6 updates. 
For A = 0.45, R(t) is constant after 2 x 10 5 updates, 
while the domain size for A — 0.25 is not constant even 
after 3 x 10 6 updates. However, we expect that for even 
small values of A, R(t) becomes constant. This behavior 
of R(t) is called pinning and is a hallmark of quenched 
disorder. As we increase the amplitude of the quenched 
disorder (increase the concentration of impurities), the 
domain growth stops and eventually becomes pinned. As 
Fig. ||b shows, in the case of T > the domains become 
larger for each value of A, but no pinning behavior is ob- 
served. In fact, the domains continue to grow throughout 
the duration of the simulation, as expected from ||. In 
both cases, two distinct time regimes are found: a rapid 
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growth in the domain size during the first 10 3 — 10 4 up- 
dates, followed by a slower growth after some crossover 
region, characterized by a time t x , that is dependent on 
the amplitude of the quenched disor der. 

Since the scaling hypothesis, Eq. (2.14), was obtained 
from considering the late stage growth behavior where 
t ^S> t x , we are interested in the scaling behavior of the 
asymptotic domain size. Fig. ^ shows the asymptotic val- 
ues Rf for each value of A for (a) T = and (b) T > 0. 
We find a linear relation, confirming the first part of the 
scaling hypothesis, namely that the domain size at late 
times scales as R ~ A~P with /3 ~ 1. We also consider 
scaling of the crossover time, t x , for each A. We assume 
that the c rosso ver time scales as t x ~ A~ 7 , where 7 ~ 3 
from Eq. (|2.15|) . Fig. shows the simulation results for 
R(t) rescaled by Eq. (§jj) with = 1,7 = 3. 

Notice that the zero temperature result, shown in 
Fig. |^a, collapses quite well over four decades in rescaled 
time, clearly confirming the scaling hypothesis. The fi- 
nite temperature result, Fig. |^b, scales well over the 
early time region only, with a clear breakdown in scal- 
ing at la te times. This suggests that the exponent 9 in 



Eq. (1.1) has a functional dependence on the amplitude 
of the quenched disorder and the temperature 

= 0(A,T). (3.1) 

We discuss reasons for this in the following section. 



IV. DISCUSSION 

We find that a time-dependent Ginzburg-Landau 
model with quenched disorder is capable of reproduc- 
ing similar results to those found in earlier Monte Carlo 
simulations ||-^) , more recent cell dynamical simulations 
and other TDGL simulations (T||. Our main in- 
terest here is the determination of the scaling function 
relating the domain size R to the amplitude of the disor- 
der A. For T = 0, we find scaling similar to that proposed 
by Grest and Srolovitz [0 as well as Glotzer et. al. ||] 
in MC simulations. Grest and Srolovitz find a = 1/2 
while Glotzer et. al. find a = 1/4 rather than the 1/3 
observed in our TDGL simulations. The reason for the 
different values of a is that MC simulations rarely see the 
predicted Lifshitz-Slyozov growth exponent of 1/3 since 
particle diffusion in an MC simulation takes many more 
updates than in a TDGL simulation. We expect that at 
low temperatures a MC simulation using an accelerated 
algorithm 20 1 or simply allowed to run sufficiently long 
will find a similar scaling function with a — 1/3. 

Nevertheless, comparison of the scaling functions for 
the zero temperature models described above enables us 
to find a relationship between the phenomenological dis- 
order parameter A and the concentration c of quenched 
impurities used in Monte Carlo simulations. The scaling 



simulations is given by Eq. (2.14) as R — A~P f(tA~<) 
with /3 ~ 1. Since we expect the domain size to scale the 
same way in both models, we anticipate the relationship 
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(4.1) 



between the amplitude of the quenched disorder and the 
impurity concentration. This result indicates that direct 
MC simulations can be replaced by TDGL simulations 
or even more efficient cell dynamical simulations, since 
the relationship between the amplitude of the quenched 
disorder can now be understood directly in terms of the 
impurity concentration for a given model. 

We note that the scaling found in the simulations dis- 
cussed above was observed either at zero temperature || 
or effectively at zero temperature |J. As our results in 
Fig. (JHJd) show, the scaling behavior for quenched dis- 
order models at finite temperature is more complicated. 
Puri and Parekh jll| observe that the asymptotic behav- 
ior of their cell dynamical model was not independent 
of the disorder parameter A. Similarly, the MC results 
of Grest and Srolovitz at T > also suggests that the 
asymptotic behavior depends on the concentration of im- 
purities. Thus we do not expect that the simple scaling 
hypothesis, Eq. (2.14), would be valid for T > 0. For late 
stag e do main growth, we expect that the exponent 8, in 
Eq. (1.1), is a function of A and T. However, a scaling 
function describing both the finite and zero temperature 
behavior of quenched disorder models might still exist. 
A detailed study of one or more of these models (MC, 
TDGL, or cell dynamical) spanning a range of tempera- 
ture and disorder parameters would likely shed consider- 
able light on this matter. Such a study would involve a 
large computational effort, and is well beyond the scope 
of the present work. 

Finally, we suggest that the type of scaling behavior 
found in a variety of quenched disorder models should 
be experimentally observable. Indeed, experiments could 
test the validity of the TDGL model discussed here, and 
could provide interesting insights into the universality of 
the scaling we report. 

We thank L. A. N. Amaral, A. Coniglio, R. Cuerno, 
S. C. Glotzer, S. Havlin, S. Sastry and G. Valley for 
helpful discussions. The Center for Polymer Studies is 
supported in part by funds from the NSF, and STH is 
supported by an NIH doctoral training fellowship. All 
simulations were done using the 32-node CM-5 at Boston 
University's Center for Computational Science. 



form suggested by MC simulations is R 



,-0' 



f(tC* ) 



with f3' = 1/d tel. The scaling form suggested by our 
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FIG. 1. Typical simulation results showing the domains for 
various values of A after 3 x 10 6 updates. The top row is the 
case T = and the bottom row is T > 0. Note the decreas- 
ing domain size as A increases. Also, the interface between 
the two phases is thicker and rougher for T > compared to 
T = 0. 

FIG. 2. (a) Spherically-averaged structure factor S(k, t) vs. 
wave number k for A — at various times. The peak value 
increases with time and approaches k = 0. (b) S(k, t) vs. k 
for A = 0.45. Note that S(k, t) is smaller, and k max is larger, 
than in (a). 

FIG. 3. (a) Double logarithmic plot of the average domain 
size R(t) vs. time t at T = for several values of A. The 
pinning regime is clearly reached for all disorder amplitudes 
except, possibly, for A = 0.25. The t 1 ^ 3 behavior is apparent 
at short times, while R ~ A^ 1 is found in the late stages, (b) 
R(t) vs. t for T > for several values of A. Note that pinning 
is not observed for any A . 

FIG. 4. (a) Inverse of the asymptotic domain size 1/ Rf vs. 
disorder amplitude A for T — 0. Rf is taken to be the average 
domain size after 3 x 10 6 updates. The linear relation con- 
firms the scaling hypothesis Eq. (2.13) with (3 — 1. (b) 1/Rf 
vs. A for T > 0. Note the deviations from a linear relation. 

FIG. 5. (a) Scaled domain size RA vs. scaled time tA 3 
for T — 0, using f) = 1 and 7 = 3 in the scaling relation 
Eq. (2.13). The collapse onto a single curve further confirms 
our scaling hypothesis, (b) Same plot as (a) using simulations 
with T > 0. 



TABLE I. Asymptotic values, Rf, of the domain size for 
various values of the quenched disorder as calculated from the 
simulations for both T = and T > 0. We used a 256x256 
lattice with spacing Ax = 0.30. The evolution of Eq. (2.9) 
was performed for 3xl0 6 updates with a time interval of 0.15. 
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